##################################################################################################
#Ms: Emergency Powers in Response to COVID-19: Policy diffusion, Democracy, and Preparedness
#Authors: Lundgren et al
#Date: September 2020
#Software: R 3.6.3
#Hardware: Mac
#Data available via Harvard Dataverse
##################################################################################################


# Read required libraries 
# Use "install.libraries" function as required

library(stargazer)
library(lme4)
library(survival)
library(rms)

# Read in data (set source directory to where data files are saved)

load("~/Dropbox/Covid paper/Data/data_cs.RData")
load("~/Dropbox/Covid paper/Data/data_tscs.RData")


# Figure 2

temp<-data_tscs[data_tscs$soe_s==1,]
temp<-temp[!is.na(temp$iso_code),]

my.surv.object <- Surv (temp$time,temp$soe_s)
fit<-npsurv(my.surv.object~1)
par(mfrow=c(1,1))
survplot(fit,xlab="Days since January 1, 2020",ylab="Proportion of states with SOE",lwd=3,conf="none",xlim=c(0,150),fun=function(y) 1-y,ylim=c(0,1))
text(0,.52, "31 December 2019: China reports pneumonia cases in Wunan to WHO",srt=90,cex=.6)
segments(0,0,0,.20)
text(71,.75, "11 March 2020: WHO declares COVID-19 as a pandemic",srt=90,cex=.6)
segments(71,0.04,71,.49)
text(30,.58, "30 January 2020: WHO declares COVID-19 a PHEIC under the IHR",srt=90,cex=.6)
segments(30,0,30,.27)


# Figure 3

mod1<-glm(soe~soe_prop+vdem_libdem+I(vdem_libdem^2)+durable+ghs_index+historical_soe_binary_imput+gdp_per_capita+who_let+wdi_pop65+population_density,data=data_cs,family="binomial")

par(mfrow=c(1,3))

newdata1 <- with(data_cs, data.frame(soe_prop=seq(0,.8,.01),vdem_libdem=mean(na.omit(vdem_libdem)), durable=mean(na.omit(durable)),ghs_index=mean(na.omit(ghs_index)),historical_soe_binary_imput=0,gdp_per_capita=mean(na.omit(gdp_per_capita)),who_let=mean(na.omit(who_let)),wdi_pop65=mean(na.omit(wdi_pop65)),population_density=mean(na.omit(population_density))))
newdata1$rankP <- predict(mod1, newdata = newdata1, type = "response")
newdata2 <- cbind(newdata1, predict(mod1, newdata = newdata1, type = "link",se = TRUE))
newdata3 <- within(newdata2, {PredictedProb <- plogis(fit)
LL <- plogis(fit - (1.96 * se.fit))
UL <- plogis(fit + (1.96 * se.fit))
})
plot(newdata3$soe_prop,newdata3$PredictedProb,type="l",lwd=2,ylim=c(0,1),ylab="Predicted probability of SOE",xlab="Proportion of states in region with SOE")
lines(newdata3$soe_prop,newdata3$UL,type="l",lwd=2,lty=3)
lines(newdata3$soe_prop,newdata3$LL,type="l",lwd=2,lty=3)
legend("topleft",c("Predicted probability","95% CI"),lty=c(1,3),lwd=2,bty="n")

newdata1 <- with(data_cs, data.frame(soe_prop=mean(na.omit(data_cs$soe_prop)),vdem_libdem=seq(0,.9,.01), durable=mean(na.omit(durable)),ghs_index=mean(na.omit(ghs_index)),historical_soe_binary_imput=0,gdp_per_capita=mean(na.omit(gdp_per_capita)),who_let=mean(na.omit(who_let)),wdi_pop65=mean(na.omit(wdi_pop65)),population_density=mean(na.omit(population_density))))
newdata1$rankP <- predict(mod1, newdata = newdata1, type = "response")
newdata2 <- cbind(newdata1, predict(mod1, newdata = newdata1, type = "link",se = TRUE))
newdata3 <- within(newdata2, {PredictedProb <- plogis(fit)
LL <- plogis(fit - (1.96 * se.fit))
UL <- plogis(fit + (1.96 * se.fit))
})
plot(newdata3$vdem_libdem,newdata3$PredictedProb,type="l",lwd=2,ylim=c(0,1),ylab="Predicted probability of SOE",xlab="Liberal democracy index")
lines(newdata3$vdem_libdem,newdata3$UL,type="l",lwd=2,lty=3)
lines(newdata3$vdem_libdem,newdata3$LL,type="l",lwd=2,lty=3)
legend("topleft",c("Predicted probability","95% CI"),lty=c(1,3),lwd=2,bty="n")

newdata1 <- with(data_cs, data.frame(soe_prop=mean(na.omit(soe_prop)),vdem_libdem=mean(na.omit(vdem_libdem)), durable=mean(na.omit(durable)),ghs_index=seq(0,85,1),historical_soe_binary_imput=0,gdp_per_capita=mean(na.omit(gdp_per_capita)),who_let=mean(na.omit(who_let)),wdi_pop65=mean(na.omit(wdi_pop65)),population_density=mean(na.omit(population_density))))
newdata1$rankP <- predict(mod1, newdata = newdata1, type = "response")
newdata2 <- cbind(newdata1, predict(mod1, newdata = newdata1, type = "link",se = TRUE))
newdata3 <- within(newdata2, {PredictedProb <- plogis(fit)
LL <- plogis(fit - (1.96 * se.fit))
UL <- plogis(fit + (1.96 * se.fit))
})
plot(newdata3$ghs_index,newdata3$PredictedProb,type="l",lwd=2,ylim=c(0,1),ylab="Predicted probability of SOE",xlab="Health security index")
lines(newdata3$ghs_index,newdata3$UL,type="l",lwd=2,lty=3)
lines(newdata3$ghs_index,newdata3$LL,type="l",lwd=2,lty=3)
legend("topleft",c("Predicted probability","95% CI"),lty=c(1,3),lwd=2,bty="n")



# Figure 4


temp<-aggregate(total_deaths~date+region,data=data_tscs,FUN=sum)

par(mfrow=c(2,4))
for (i in unique(temp$region[temp$region!="South Asia"])){
  temp2<-temp[temp$region==i,]
  plot(seq(as.Date("2020-01-01"),as.Date("2020-06-12"),1),log(temp2$total_deaths),main=i,type="l",lwd=1.5,ylim=c(0,12.5),xlab="Date",ylab="Total deaths in region (log)") 
  temp3<-aggregate(soe_s~date,data=data_tscs[data_tscs$soe_s==1 & data_tscs$region==i,],FUN=sum)
  for (k in 1:length(temp3$date)){
    temp3$deaths[k]<-temp2$total_deaths[temp2$date==temp3$date[k]]
    print(k)
  }
  points(as.Date(temp3$date),log(temp3$deaths),pch=1,cex=temp3$soe_s)
}
temp2<-temp[temp$region=="South Asia",]
plot(seq(as.Date("2020-01-01"),as.Date("2020-06-12"),1),log(temp2$total_deaths),main="South Asia",type="l",lwd=1.5,ylim=c(0,12.5),xlab="Date",ylab="Total deaths in region (log)") 



# ONLINE APPENDIX

# Table A1 Descriptive statistics, cross-sectional data structure

datasummary<-subset(data_cs, select=c("soe","soe_prop","vdem_libdem","durable","ghs_index","historical_soe_binary_imput","gdp_per_capita","who_let","wdi_pop65","population_density"))
stargazer(datasummary, digits=2,type="text")


# Table A2 Descriptive statistics, tscs data structure

datasummary<-subset(data_tscs, select=c("soe_s","lag.cum_soe_region","lag.cum_soe","vdem_libdem","polity2","durable","ghs_index","lag2.StringencyIndex","wdi_pop65","wdi_mortinf","gdp_per_capita","population_density","total_deaths_per_million","em","historical_soe_binary_imput"))
stargazer(datasummary, digits=2,type="text")



# Table A3 TSCS models 
pvars <- c("lag.cum_soe","vdem_libdem",
           "durable","ghs_index",
           "wdi_pop65","wdi_mortinf","gdp_per_capita","population_density","total_deaths_per_million","time","lag.StringencyIndex","lag2.StringencyIndex","historical_soe_binary_imput","lag.cum_soe_region")
datsc <- data_tscs
datsc[pvars] <- lapply(datsc[pvars],scale)

# Model 1: main model
model1 <- glmer(soe_s ~  lag.cum_soe_region+vdem_libdem+I(vdem_libdem^2)+log(durable+1)+ghs_index+lag2.StringencyIndex+wdi_pop65+wdi_mortinf+gdp_per_capita+population_density+total_deaths_per_million+time+ I(time^2)+I(time^3)+(1  | iso_code),family=binomial(link='logit'),data=datsc,control = glmerControl(optimizer = "bobyqa",calc.derivs = F,optCtrl=list(maxfun=1e5)))

# Model 2: with constitutional em powers 
model2 <- glmer(soe_s ~  lag.cum_soe_region+vdem_libdem+I(vdem_libdem^2)+log(durable+1)+ghs_index+lag2.StringencyIndex+wdi_pop65+wdi_mortinf+gdp_per_capita+population_density+total_deaths_per_million+em+time+ I(time^2)+I(time^3)+(1  | iso_code),family=binomial(link='logit'),data=datsc,control = glmerControl(optimizer = "bobyqa",calc.derivs = F,optCtrl=list(maxfun=1e5)))

# Model 3: with historical soes 
model3 <- glmer(soe_s ~  lag.cum_soe_region+vdem_libdem+I(vdem_libdem^2)+log(durable+1)+ghs_index+lag2.StringencyIndex+wdi_pop65+wdi_mortinf+gdp_per_capita+population_density+total_deaths_per_million+historical_soe_binary_imput+time+ I(time^2)+I(time^3)+(1  | iso_code),family=binomial(link='logit'),data=datsc,control = glmerControl(optimizer = "bobyqa",calc.derivs = F,optCtrl=list(maxfun=1e5)))

# Model 4: with polity data instead of v-dem
model4 <- glmer(soe_s ~  lag.cum_soe_region+polity2+I(polity2^2)+log(durable+1)+ghs_index+lag2.StringencyIndex+wdi_pop65+wdi_mortinf+gdp_per_capita+population_density+total_deaths_per_million+time+ I(time^2)+I(time^3)+(1  | iso_code),family=binomial(link='logit'),data=datsc,control = glmerControl(optimizer = "bobyqa",calc.derivs = F,optCtrl=list(maxfun=1e5)))

# Model 5: without the united states 
datsc2<-datsc[datsc$iso_code!="USA",]
model5 <- glmer(soe_s ~  lag.cum_soe_region+vdem_libdem+I(vdem_libdem^2)+log(durable+1)+ghs_index+lag2.StringencyIndex+wdi_pop65+wdi_mortinf+gdp_per_capita+population_density+total_deaths_per_million++time+ I(time^2)+I(time^3)+(1  | iso_code),family=binomial(link='logit'),data=datsc2,control = glmerControl(optimizer = "bobyqa",calc.derivs = F,optCtrl=list(maxfun=1e5)))

stargazer(model1,model2,model3,model4,model5,type="text",digits=2,df=F)
stargazer(model1,model2,model3,model4,model5,type="html",digits=2,df=F)

